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We compute exact asymptotic results for the probability of the occurrence of large deviations of 
the largest (smallest) eigenvalue of random matrices belonging to the Gaussian orthogonal, unitary 
and symplectic ensembles. In particular, we show that the probability that all the eigenvalues of an 
(A'^ X A'^) random matrix are positive (negative) decreases for large N as ^ exp[—Pd{0)N'^] where 
the Dyson index l3 characterizes the ensemble and the exponent 0{0) = (In3)/4 = 0.274653 ... is 
universal. We compute the probability that the eigenvalues lie in the interval [^i , Ca] which allows 
us to calculate the joint probability distribution of the minimum and the maximum eigenvalue. 
As a byproduct, we also obtain exactly the average density of states in Gaussian ensembles whose 
eigenvalues are restricted to lie in the interval [1^1,(^2], thus generalizing the celebrated Wigner semi- 
circle law to these restricted ensembles. It is found that the density of states generically exhibits 
an inverse square-root singularity at the location of the barriers. These results are confirmed by 
numerical simulations. 

PACS numbers: 02.50.-r, 02.50.Sk, 02.10.Yn, 24.60.-k, 21.10.Ft 

I. INTRODUCTION 

Studies of the statistics of the eigenvalues of random matrices have a long history going back to the seminal work 
of Wigner [ij. Random matrix theory has been successfully applied in various branches of physics and mathematics, 
including in subjects ranging from nuclear physics, quantum chaos, disordered systems, string theory and even in 
number theory [3]. Of particular importance are Gaussian random matrices whose entries are independent Gaussian 
variables Depending on the physical symmetries of the problem, three classes of matrices with Gaussian entries 
arise {N x TV) real symmetric (Gaussian Orthogonal Ensemble (GOE)), (N x N) complex Hermitian (Gaussian 
Unitary Ensemble (GUE)) and (27V x 2A^) self-dual Hermitian matrices (Gaussian Symplectic Ensemble (GSE)). In 
these models the probability distribution for a matrix M in the ensemble is given by 

KM)(xexp(^-^(M,M)^ , (1) 

where (M, M) is the inner product on the space of matrices invariant, under orthogonal, unitary and symplectic 
transformations respectively and the parameter /3 is the Dyson index. In these three cases the inner products and the 
Dyson indices are given by 

(M,M) = Tr(M2); /3 = 1 GOE (2) 

(M,M) = Tr(M*M); /3 = 2 GUE (3) 
(M,M) = Tr(Ml'M); /3 = 4 GSE (4) 

Jv>( , where •* denotes the hermitian conjugate of complex valued matrices and denotes the symplectic conjugate on 
5_j quaternion valued matrices. The above quadratic actions are the simplest forms (corresponding to free fields) of 
^ matrix models which have been extensively studied in the context of particle physics and field theory. 

A central result in the theory of random matrices is the celebrated Wigner semi-circle law. It states that for large 

N and on an average, the N eigenvalues lie within a finite interval —\/2N, \/2N , often referred to as the Wigner 
'sea'. Within this sea, the average density of states has a semi-circular form (see Fig. ([!])) that vanishes at the two 
edges -V2iV and y/2N 



2iV 



1/2 



(5) 



The above result means that, if one looks at the density of states of a typical system described by one of the three 
ensembles above, for a large enough system, it will resemble closely the Wigner semi-circle law. 
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FIG. 1: The dashed Une shows the Wigner semi-circular form of the average density of states. The largest eigenvalue is centered 
around its mean y/2N and fluctuates over a scale of width N^^^^ . The probability of fluctuations on this scale is described by 
the Tracy- Widom distribution (shown schematically). 



While the semi-circle law provides a global information about how the eigenvalues are typically distributed, unfor- 
tunately it does not contain enough information about the probabilities of rare events. The questions concerning rare 
events have recently come up in different contexts. For example, string theorists have recently been confronted by 
the possibility that there may be a huge number of effective theories describing our universe. The landscape made up 
of these theories is called the string landscape. This seemingly embarrassing situation may however help to explain 
certain fine tuning puzzles in particle physics. The basic argument is as follows. Our universe is one which supports 
intelligent life and this requires that the ratios of certain fundamental constants lie in specific ranges. The fine tuning 
we observe is thus a necessary condition that we are there to describe it. Other possible universes would have different 
vacua but there would be no intelligent life to study them. This approach to string theory is called the anthropic prin- 
ciple based string theory and is a subject of current and intense debate. In 0] the authors carried out an analysis 
of the string landscape based solely on the basis that it is described by a large N multi-component scalar potential. 
Of particular interest is the determination of the typical properties of the vacua in the string landscape based only on 
assumptions about the dimensionality of the landscape and other simple general features. The motivation for these 
studies is to determine to what extent the string landscape is determined by large N statistics and what features 
depend on the actual structure of the underlying string theory. 

One of the main questions posed in Q is: for a {N x N) Gaussian random matrix, what is the probability Pn that 
all its eigenvalues are positive (or negative) 

Pjv = Prob[Ai > 0,A2 > 0, . ..,Aiv > 0]? (6) 

From the semi-circle law, we know that on an average half the eigenvalues are positive and half of them are negative. 
Thus the event that all eigenvalues are negative (or positive) is clearly an atypical rare event. This question arises in 
the so called 'counting problem' of local minima in a random multi-field potential or a landscape. Given a stationary 
point of the landscape, if all eigenvalues of the Hessian matrix of the potential are positive, clearly the stationary 
point is a local minimum. Thus, the probability that all eigenvalues of a random Hessian matrix are positive provides 
as estimate for the fraction of local minima amongst the stationary points of the landscape. In particular, the authors 
of [3| studied the case where the Hessian matrix was drawn from a GOE ensemble (/3 = 1). Thus, in this context Pjv 
is just the probability that a random GOE matrix is positive definite. This probability has also been studied in the 
mathematics literature Q and one can easily compute Pjv for smaller values of = 1,2,3. For example, one can 
show that [5i] 

The interesting question is how Pn behaves for large A^? It was argued in Q that for large iV, Pn decays as 
P/v ^ exp (— 6'(0)iV^) where the decay constant 0(0) was estimated to be « 1/4 numerically and via a heuristic 
argument. The scaling of this probability with N"^ is not surprising and has been alluded to in the literature for a 
number of years, notably in relation to studies of the distribution of the index (the number of negative eigenvalues) of 
Gaussian matrices @, 0|- However an exact expression for 0(0) was not available until only recently, when the short 
form of this paper [8| was published. In Q we had shown that for all the three Gaussian ensembles, to leading order 
in large N 

1 3 

P/v ~ exp[-/36'(0)A^2j. .^^^ere 0(0) = ^ = (0.274653 . . .). (8) 
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Interestingly, the probability Pj^ has also recently shown up in a rather different problem in mathematics. Dedieu and 
Maljovich has shown [5] recently that Ppf is exactly equal to the expected number of minima of a random polynomial 
of degree at most 2 and N variables. Our result in ^ thus provided an exact answer to this problem for large N 

To put our results in a more general context, we note that the semi-circle law tells us that the average of the maximum 
(minimum) eigenvalue is \/2N [-\/2N). However, for finite but large N, the maximum eigenvalue fluctuates, around 
its mean \/2N, from one sample to another. Relatively recently Tracy and Widom Q proved that these fluctuations 
typically occur over a narrow scale of ~ 0(iV~^/^) around the upper edge \/2N of the Wigner sea (see Fig. 1). 
More precisely, they showed that asymptotically for large N, the scaling variable ^ = -^2-^^^^ [Amax — V2N] has 
a limiting A^-independent probability distribution, Prob[^ < x] — Ffj{x) whose form depends on the value of the 
parameter /3 = 1, 2 and 4 characterizing respectively the GOE, GUE and GSE. The function Fp{x), computed as a 
solution of a nonlinear differential equation |9|], approaches to 1 as x — > oo and decays rapidly to zero as x — *■ — oo. 
For example, for P — 2, F2{x) has the following tails Q, 

F2{x) 1 - O (^exp[-4a;^/V3]) as a; ^ oo 

exp[-|a;|^/12] as x ^ -oo. (9) 

The probability density function dFfj / dx thus has highly asymmetric tails. The distribution of the minimum eigenvalue 
simply follows from the fact that Prob[Amin ^ C] = Prob[Amax < — CJ - Amazingly, the Tracy- Widom distribution 
has since emerged in a number of seemingly unrelated problems [10| such as the longest increasin g su bsequence 
problem [llj, directed polymers in (1 -I- l)-dimensions various (1-1- l)-dimensional growth models |l3l| . a class of 
sequence alignment problems [l^ . mesocopic fluctuations in dity metal grains and semiconductor quantum dots [l5l | 
and also in finance 'l6]. 

The Tracy- Widom distribution describes the probability of typical and small fluctuations of Amax over a very narrow 
region of width ^ 0{N~'^/^) around the mean (A,„ax) ~ V2N. A natural question is how to describe the probability 
of atypical and large fluctuations of Xmax around its mean, say over a wider region of width ~ 0(7Vi/2)? For example, 
the probability Pjv that all eigenvalues are negative (or positive) is the same as the probability that Amax < 
(or Amin > 0). Since (Amax) ~ \/2N, this requires the computation of the probability of an extremely rare event 
characterizing a large deviation of ~ —0{N^/^) to the left of the mean. In we calculated the exact large deviation 
function associated with large fluctuations of ~ —0{N^^'^) of Amax to the left of its mean value ^/2N. It was shown 
that for large N and for all ensembles 



Prob [Amax <t,N]r^ exp 



/2N~t 



N 



(10) 



where t ^ 0{N^/^) < \/2N is located deep inside the Wigner sea. The large deviation function <&(?/) is zero for y < 0, 
but is nontrivial for y > which was computed exactly in Q . For small deviations to the left of the mean, taking the 
?/ — !■ limit of ^{y), one recovers the left tail of the Tracy- Widom distribution as in Eq. ([5]). Thus our result for large 
deviations of ^ —0{N'^/'^) to the left of the mean is complementary to the Tracy- Widom result for smo^/ fluctuations 
of ^ — 0(iV^^/^) and the two solutions match smoothly. Also, the probability P/v that all eigenvalues are negative 
(or positive) simply follows from the general result in Eq. (jlOp by putting t ^ Q, 

Pn - Prob [Amax < 0, TV] ~ exp [Vt] N^] (11) 



thus identifying 9{Q) = <I> (V2) — (ln3)/4, a special case of the general large deviation function. 

The purpose of this paper is to provide a detailed derivation of the above results announced in S], as well as 
numerical results in support of our analytical formulas. In addition, we also derive asymptotic results for the joint 
probability distribution of the minimum and the maximum eigenvalue. 

Statistical analysis motivated by anthropic considerations in string theory or random polynomials in mathematics 
may seem a long way from laboratory based physics, however similar questions appear naturally also in providing 
criteria of physical stability in dynamical systems or ecosystems (TtI . [isj . Near a fixed point of a dynamical system, 
one can linearize the equations of motion and the eigenvalues of the corresponding matrix associated with the linear 
equations provide important informations about the stability of the fixed point. For example, if all the eigenvalues 
are negative (or positive) the fixed point is a stable (or unstable) one. In this context, another important question 
arises naturally. Suppose that the dynamical system is close to stable (unstable) fixed point, i.e., all the eigenvalues 
are negative (positive). Given this fact, one may further want to know how these negative (positive) eigenvalues are 
distributed. In other words, what is the average density of states of the negative (or positive) eigenvalues given the 
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fact that one is close to a stable (or unstable) fixed point. In this paper we will calculate the density of states in this 
conditioned ensemble and we will see that it is quite different to the Wigner semi-circle law. 

Recently the problem of determining the stability of the critical points of Gaussian random fields in large dimensional 
spaces was analyzed [l^ [l^l . The Hessian matrix in this case does not have the statistics of a Gaussian ensemble and 
the probability that a randomly chosen critical point is a minimum case be shown to decay as exp^—Nip) and the 
exponent ip can be explicitly calculated in terms of the two point correlation function of the Gaussian field. Thus the 
scaling in N is quite different to the random matrix case and this scaling obviously makes minima much more likely 
and yields a more usual thermodynamic scaling of the entropy of critical points 19]. The statistics of the Gaussian 
field problem are perhaps more relevant to statistical landscape scenarios in string theory. 

The paper is organized as follows. In Section II we begin by recalling the Goulomb gas representation of the 
distribution of eigenvalues of Gaussian matrices and show how our problem can be formulated by placing a hard wall 
constraint on the Coulomb gas. This is the key step in the method and the technique has since been applied to analyze 
the statistics of critical points of Gaussian random fields [l^ and also to study the probability of rare fluctuations of 
the maximal eigenvalue of Wishart random matrices [llj . We then show how in the large N limit the problem can be 
solved using a saddle point computation of a functional integral and wc discuss the features of our analytic results. 
In Section III, we extend this method to compute the asymptotic joint probability distribution of the minimum and 
the maximum eigenvalue. This requires studying the Coulomb gas confined between two hard walls. In Section IV we 
carry out some numerical work to confirm our predictions about the probability of extreme deviations of the maximal 
eigenvalue. We show how the Coulomb gas formulation can again be exploited even for numerical purposes. Finally 
we present our conclusions in Section V. 



II. THE COULOMB GAS FORMULATION AND THE PROBABILITY OF RARE FLUCTUATIONS 

The joint probability density function (pdf) of the eigenvalues of an x iV Gaussian matrix is given by the classic 
result of Wigner [ll, 



P(Ai, A2, . . . , Aat) = Sat exp 



N 

^A?-^ln(|A,-A,|) 

^1 *7^J 



(12) 



where Bn normalizes the pdf and f3 — 1, 2 and 4 correspond respectively to the GOE, CUE and GSE. The joint law 
in Eq. (fT2|) allows one to interpret the eigenvalues as the positions of charged particles, repelling each other via a 2-d 
Coulomb potential (logarithmic); they are confined on a 1-d line and each is subject to an external harmonic potential. 
The parameter f3 that characterizes the type of ensemble can then be interpreted as the inverse temperature. 

Once the joint pdf is known explicitly, other statistical properties of a random matrix can, in principle, be derived 
from this joint pdf. In practice, however this is often a technically daunting task. For example, suppose we want to 
compute the average density of states of the eigenvalues defined as p{X,N) — J2iLi{'^i^ ^ ^i))/^^ which counts the 
average number of eigenvalues between A and A + c?A per unit length. The angled bracket () denotes an average over 
the joint pdf. It then follows that p(A, N) is simply the marginal of the joint pdf, i.e, we fix one of the eigenvalues 
(say the first one) at A and integrate the joint pdf over the rest of the {N — 1) variables. 

. N .00 N 

p{\ N)^- Y,{5{\ - K)) = / H dX, P(A, A2, . . . , Aat). (13) 

i=l "'-°Oi=2 

Wigner computed this marginal and showed [l[ that for large N and for all /3 it has the semi-circular form in Eq. (O . 

Here we are interested in calculating the probability QNiC) that all eigenvalues are greater than some value C- This 
probability is clearly also equal to the cumulative probability that the minimum eigenvalue Amin = min(Ai, A2, . . . , Ajv) 
is greater than i.e., 

Qjv(C) = Prob [Ai > C, A2 > C, • • • , Aat > C] = Prob[A,„i„ > C, (14) 

Since the Gaussian random matrix has the x —x symmetry, it follows that the maximum eigenvalue Amax has the 
same statistics as — Amin- Hence, it follows that 

Prob[A,„ax <t,N]^ Qn{C = -t). (15) 

Hence, knowing QNiC) '^iU f^lso allow us to compute the cumulative distribution of the maximum. Also, note that 
the probability Pjv that all eigenvalues are positive (or negative), as defined in the introduction, is simply 

Pn = QNiO). (16) 
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In what follows, we will compute QNiC) in the scaling limit where C '-^ yiV for large N. For this we will employ the 
saddle point method in the framework of the Coulomb gas. 
By definition 

/OO /"OO 
...J^ P(Ai,A2,...,Aw)dAidA2...dA^ (17) 

where P is the joint pdf in Eq. (|12p. This multiple integral can be written as a ratio 

QNiO - /f^ (18) 

Zjv(-oo) 

where the partition function Zn{C,) is defined as 

exp -- I 



Zm{0 = / . . . / exp 



^A,?-^ln(|A.-A, 



dXidX2 . . .dXN- (19) 



Note that the normalization constant = l/Zff{—oo). Clearly, for any finite C, Zi^{C.) represents the partition 
function of the Coulomb gas which is constrained in the region [C, oo], i.e., it has a hard wall at ^ ensuring that there 
are no eigenvalues to the left of 

From the semi-circle law, it is evident a typical eigenvalue scales as ^/N for large A''. This is also evident from Eq. 
where the Coulomb interaction term typically scales as N'^ for large TV while the energy corresponding to the 
external potential scales as A^A^. In order that they balance, it follows that typically A ^ ^/N. It is therefore natural 
to rescale /i^ — Xi/y/N so that the rescaled eigenvalue fii ^ 0(1) for large N. In terms of the rescaled eigenvalues, 
the partition function reads 

TV 



N 



Zn{0'x / Y[dfi,exp{-l3H{fi)) 

Jui>-i^ ._i 

where the Hamiltonian of the Coulomb gas is 

= y Ea*' - ^ J2H\^'^ - M,l). (20) 

i=l i^j 

A. EYinctional Integral, Large N Saddle Point Analysis and the Constrained Charge Density 

We first define the normalized (to unity) spatial density field of the particles fii as 

1 ^ 

pi^')^JJY.^(^'-^'^y (21) 



This is just a 'counting' function so that p{fi) dfi counts the fraction of eigenvalues between /x and fi + dfi. The energy 
of a configuration of the fii 's can be expressed in terms of the density p as 

H[p] = N^£[p] (22) 

where 

^ip] = l J dpLp^-" p{p)~]^ J dpd^i'pip)p{^i')H\fI- fi'\) + ^ J dp p{p)\n{l{p)) . (23) 

The last term above removes the self interaction energy from the penultimate term and l{p) represents a position 
dependent cut-off. Dyson argued that l{p) ^ l/p{p) and hence this correction term has the form 

dpi p{p) In {l{p)) ^ - I dp p{p) In (p(/x)) -I- C, (24) 
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where C" is a constant that cannot be determined by this argument, but one may assume that it is independent of 
the position of the hard wall. However, in this paper, we will be interested only in the leading 0{N'^) behavior and 
hence precise value of the constant C does not matter. 

The hard wall constraint can now be implemented simply by the condition 

pU) = for /i < (25) 

The partition function may now be written as a functional integral over the density field p as 

Zn{0 = j d[p] J[p] exp {-PN^£[p]) , (26) 

where J[p\ is the Jacobian involved in changing from the coordinates p,i to the density field p. Physically this Jacobian 
takes into account the entropy associated with the density field p. For the sake of completeness we will re-derive a 
familiar form of the Jacobian J[p]. Clearly J can be written, up to a constant prefactor D^r, as 



J[p\=Dn l\{dp,5 
•' i=i 



(27) 



where the integration range of the pi above are restricted to the appropriate region. One now proceeds by making 
a functional Fourier transform representation of the delta function at each point p. This gives 



J[p\ ^ D'j^ JYI '^t^i'^id] exp 



dp g{p) 



Np{p)-Y,5{p-pi) 



(28) 



where each g{p) integral is along the imaginary axis and D'^ is a constant prefactor. The integral over the pi may 
now be carried out giving 



J[p\ = D'n j d[g\ exp N j dp g[p) p{p) J]^ d/i^ exp [-^(/i,)] 



D'n j d[g] exp 



N / dp g{p)p{p) + N In { / exp(-.g(^)) 



(29) 



The above functional integral over g can be evaluated by saddle point for large N and the corresponding saddle point 
equation (obtained via stationarity with respect to g) is 



exp {-g{p)) 
Jdp'exp^-gip'))' 



We see that the normalization 



dp pip) = 1 

is respected. Substituting in this solution for g we find that 



N / dpp{p) In {pip}) 



J[p] = D'n exp 



from which we can read off the entropy corresponding to the density field p. 
Putting this all together we find 



(30) 



(31) 



(32) 



Zn{C) = ^a^/ ^ip]^ (J dp p{p) - l) exp ~N j dp p{p)\n{p{p)) - l3N^£[p\ 



(33) 

where is a prefactor and the delta function enforces the normalization condition in Eq. (|3ip . One simple way to 
incorporate this delta function constraint is to introduce a Lagrange multiplier C (corresponding to writing the delta 
function constraint in the Fourier representation) and rewrite the partition function as 



Zjv(C) -^W y"rfCd[p]exp 



-PN^m - N{1 -I:) I dp Pip) In ipip)) 



(34) 
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where A'j^ is a constant and the 'renormahzed' action is expHcitly 

= ^ 1 dfi fi'^ pin) - ^ 1 dndp.'p{^) p{n') \n{\fi- fi'\) + C 



dpp[p] 



(35) 



The 0{N) term in Eq. (1341) includes both the entropy of the density field and the self energy subtraction following 
the Dyson prescription discussed before. While this prescription is difficult to justify rigorously, we need not pursue 
this issue here since we are interested only in the leading O(iV^) behavior. 

The functional integral in Eq. p4|) can thus be evaluated by the saddle point method in the variable N"^ and the 
term of order N coming from the entropy and self energy subtraction is negligible in the large N limit. We thus find 
that in the saddle point analysis, for large N, 

ZN{0^exp[-PN^i:[p,]+OiN)] (36) 

where pdp) is the density field that minimizes the action S][/9] in Eq. (|35p . Minimizing it follows that Pc{p) 

satisfies the integral equation 

4 + r dp' p4p')\n{\fi^p'\), (37) 



where we have introduced the scaled variable z ~ C,/vN. The normalization and boundary conditions for Pc{p), hi 
terms of the scaled variable z, are 

/>oo 

/ dppc{p) = l; pc{p) ^ for p < z. (38) 

J z 

The partition function then reads 

Zn{z) ^eyip{-f3N'^ S{z) + 0{N)) . (39) 

where 

^(z) -min{S][p]}-E[pJ (40) 
p 

Let us remark here that in the unconstrained case (z = — oo) the Wigner semi-circle law can also be obtained from 
the saddle-point method 0- However to our knowledge the constrained problem has never been analyzed using this 
approach. In addition, the effective free energy is not extensive: it scales as N'^ rather than N due to the long-range 
nature of the logarithmic inter-particle interaction. The entropy term is extensive and is contained in the 0{N) term 
in Eq. (j39p . However the entropy is subdominant, and the free energy is dominated by the energetic component. 

Differentiating Eq. ((37)) with respect to p we get 

p=v r dp' BAin.^ (41) 

Jz M - M 

where V indicates the Cauchy principle part. It is convenient to shift the variable p by writing 

p = z + x (42) 

and where a; > is now positive and a; = denotes the location of the infinite barrier. In terms of the shifted variable, 
we denote the density field as 

Pc{p = X + z) ^ f{x;z). (43) 

Eq. (|4T|) then becomes 

x + z^V r dx' (44) 

Jo X — X 

This integral equation is of the general form 

g{x) = {H+f){x) =V r dx' (45) 
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with g{x) = X + z. The right hand side (rhs) of the above equation is just the semi-infinite Hilbert transform of the 
function f{x). The main technical challenge is to invert this half Hilbert transform. Fortunately this inversion can 
be done for arbitrary g{x) using Tricomi's theorem [2^. To apply Tricomi's theorem, we first assume, to be verified 
aposteriori, that f{x) has a finite support over [0, L], so that the integral in Eq. can be cut-off at L at the upper 
edge. Then the solution for f{x) is given by Tricomi's theorem [2^ 

fi-) = - 2 A^^ l^T'^"' I (46) 

T:-'y'x[L~x) Jo x-x J 

where C is an arbitrary constant which can be fixed by demanding that f{x) vanishes at x = L. Applying this 
formula to our case with g{x) — x + z and performing the integral in Eq. (j46p using Mathematica, we get 

1 



f{x;z) = ^L{z) - X [L{z) + 2x + 2z] , a;e[0,L(z)] 

f{x; z) = 0, X <0 and x > L{z) (47) 

where we have made the z dependence of L{z) explicit. The value of L[z) is now determined from the normalization 
of /(x; z), i.e., /q^'^"* /(x; z)dx — 1 and one gets 

L{z) = \ Wz-^ + 6 - zj (48) 

which is always positive. We also see that in terms of the variable ^ = .x + z, the support of Pci^J■) is between the 
barrier location at z and an upper edge a,t n = z + L{z) = | [vz^ + 6 -I- 2z] . 

Let us make a few observations about the constrained charge density Pc(m — x + z) = f{x] z) in Eqs. (|47|) and (j48|) . 



• Physically, the charge density f{x;z) must be positive for all x including x = 0. Now, as a; ^ 0, i.e., as one 
approaches the infinite barrier from the right, it follows from Eq. ((47|l that f{x; z) diverges as a;~^/^, i.e., the charges 
accumulate near the barrier. However, for it to remain positive at a; = 0, it follows from Eq. ((47|) that we must have 
L{z) -t- 2z > 0. This happens, using the expression of L{z) from Eq. (|48p . only for z > ~V2. Thus the result in 
Eq. ([47|) is valid only for z > — V2. To understand the significance of z > — \/2, we note that when z — —\pi, i.e., 
when the barrier is placed exactly at the leftmost edge of the Wigner sea, we recover the Wigner semi-circle law from 
Eqs. ((17)) and (^5]) . We get L{z = — \/2) = 2\/2 (the support of the full semi-circle), i.e., f{x;—^/2) is nonzero for 
< a; < 2^/2. In terms of the original variable fi — x + z = x— \/2, pdn) is nonzero in the region —V2 < p. < \f2 
and is given by, within this sea to the right of the barrier, in terms of the original variable /i 

1 



p,{p) - -V2^. (49) 

TT 

When z becomes smaller than — v^, the solution in Eq. (|47|) remains unchanged from its semi-circular form at 
z = —a/2. In other words, for z < — v^, the solution sticks to its form at z — — V2- Physically, this means that if the 
barrier is placed to the left of the lower edge of the Wigner sea, it has no effect on the charge distribution. 

• Note that the above fact, in conjunction with Eq. (|l(I)l . indicates that 

S{z) ^ S{-V2) for aU z < -V2 (50) 

including, in particular 

S{-oo) = S{-V2) (51) 

a result that we will use later. 

• The charge density pc{fJ, = x + z) = f{x;z) changes its shape in an interesting manner as one changes the 
barrier location z. For z > — \/2 the global maximum of f{x;z) is always at a; = (at the barrier) where it has a 
1/a/x integrable singularity. For — < z < —^3/4, the density is non-monotonic and in addition to the square root 
divergence at a; = 0, f{x; z) develops a local minimum and a local maximum respectively at a: = (L^V— — 8zL)/4. 
For z > — •\/3/4, the density f{x;z) decreases monotonically with increasing x. In Fig. ^ we show the behavior of 
/(a;; z) for three representative values of z. An interesting point about the eigenvalue distribution /(x; z) given Eq. 
P7)) is that, for z > — \/2 there is a strong accumulation of eigenvalues at the barrier location a; = or equivalently 
a,t n — z. In the case z = if one thinks of the eigenvalues as being associated about the vacuum of a (stable) 
field theory then there is an accumulation of modes of mass close to zero, a fact that may have consequences in the 
context of anthropic principal based string theory or in other physical systems where only stable configurations can 
be observed. 
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FIG. 2: Plots of the density of states f{x;z) as a function of the shifted variable x for z = —1 (dotted), z — (solid), and 
z = 0.5 (dashed). 

B. Large Deviations of the Maximum or the Minimum 

Having computed the constrained charged density pdfJ-), we next calculate the action S{z) = T,[pc] at the saddle 
point. In order to calculate S[pc] from Eq. we will use the explicit saddle point solution pc{n = x + z) = f{x; z) 

obtained in Eq. (|47| . In fact, one can simplify the expression of the saddle point action by using the integral 
equation ([57]) satisfied by the saddle point solution. We multiply Eq. ([57|) by Pcip) and integrate over p. Using the 
normalization J pc{p) dp = 1 one gets 



dpdp'pcip) Pc{p') ln(|M - m'I) = C 



p^ Pc{p) dp 



(52) 



where C is the Lagrange multiplier to be determined. Substituting this result in Eq. psp and the fact, J Pc{p) dp ^ 1 
we get 



1 If 

^[Pc] = --^C + - / p^pc{p)dp. 



(53) 



To determine the Lagrange multiplier C, we put p — z in Eq. ([37| . In the integral on the rhs of Eq. (|37p we then 
make the usual shift, p' = z + x' and use pdp') = f{x'; z) to get 



C 



L{z) 



dx' f{x']z)\n(x') . 



(54) 



where /(x; z) is explicitly given in Eq. (I47p and L{z) in Eq. 
Pc{p') = f{x'; z) we finally get 



Substituting the expression of C in Eq. ([55)1 and 



1 1 f^^^^ 1 
^[Pc] = ^2:^-2 J dx ln{x) f{x;z) + - J dx {x + zY f {x; z) 



L(z 



(55) 



with f{x;z) given explicitly in Eq. (I47p . The integrals can again be performed explicitly using Mathematica and 
using Eq. we get the following expression for the saddle point action 



S(z) = 

^ ' 216 



72z2 - 2z'^ 



+ (30z + 2z^) v/6 + z2 + 27 (3 + ln(1296) - 4 In (-z + ^6 + z^^ ) 



(56) 
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The partition function Zj^{(^ 



Nz) is then given by Eq. ((39|) . Note also, using Eq. (|5ip . we have 

Zn{—co) ~ exp 



(57) 



where S'(-\/2) = (3 + ln(4))/8 from Eq. Taking the ratio in Eq. ^ gives us Qjv(C), the probabihty (5iv(C) 

that all eigenvalues are to the right of ^ = -x/iV z 



'.N 



(C) = exp 



C 



0[N) 



where 9{z) = 5(2;) — S{—^/2) is given by 



1 



^ ^ 108 



36^2 - + (15z + z^)^/& + z^ + 27 (ln(18) - 2 In (-2 + ^6 + 
The probability that all eigenvalues are positive (or negative) is simply 



where 



Pn = Qiv(C = 0) « exph/3 0(0) 



0(0) = = 0.274653. 



(58) 



(59) 



(60) 



(61) 



Finally let us turn to the large deviation function associated with large negative fluctuations of ~ — 0(iV^/^) of 



Amax to the left of its mean value v2iV. Substituting the expression for Qn{0 from Eq. ([58|) in Eq. (fT5|) we get 



Prob[A„iax <t,N]^ Qn{C = -t) = exp 



0{N) 



Noting that (Amax) = V2iV, it is useful to center the distribution around the mean and rewrite Eq. ([62)1 as 



Prob[A„ 









exp 











(62) 



(63) 



where the large deviation function $(?/) — 9{y — V2) with 9{z) given explicitly in Eq. ([5^ 

One can easily work out the asymptotic behavior of $(j/) for small and large y. For example, it is easy to see that 



6V2 
2 



as y 







y 

— as y 
2 ^ 



(64) 



In particular, when v2N — t << v^, i-e, we are rather close to the right edge (Amax) — v27V of the Wigner sea (see 
Fig. 1), it follows that the scaling variable y = (\/2iV — t)/\/N « 1. Hence substituting the small y behavior of 



/Q\/2 from Eq. ([M)) in Eq. ([55]) . it follows that in this regime 



Prob[Amax <t,N]K exp 



A(|V2iVi/6(t_V27V)| 



(65) 



Note that this matches exactly with the left tail behavior of the Tracy- Widom limiting distribution for all the three 
cases /3 = 1, 2 and 4 For example, for (3 — 2, one can easily verify by comparing Eqs. (I65|) and This is to be 
expected because the Tracy- Widom distribution describes the distribution of Amax around its mean ^/2N over a scale 
~ 0(iV~^/^). If we want to investigate the probability of negative fluctuations of order {V2N — t) >> iV~^/^, we 
need to look at th e left tail of the Tracy- Widom distribution. On the other hand, those negative fluctuations of order 
N~^/^ « (■\/2N—t) << ^/N are described by the small argument behavior of the large deviation function. These two 
behaviors thus should smoothly match. As we verified above, they do indeed match smoothly, thus providing another 
confirmation of our exact result. Moreover, our large deviation function provides an alternative way to compute the 
left tail of the Tracy- Widom distribution for all /3. 
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III. COULOMB GAS BOUNDED BY TWO WALLS: THE JOINT PROBABILITY DISTRIBUTION OF 

A mill AND A max 

In this section we compute, by the Coulomb gas method, the probability Rn{CiX2) that all eigenvalues are in 
the interval [Ci,C2] where > Ci- Evidently, in the limit (2 00, i?Ar(C, 00) = Qn{C) which was computed in the 
previous section. Clearly, Rpf {(1,(2) is also the cumulative probability that Amin > Ci ^-nd Amax < C2, i.e., 

-Rjv(Ci, C2) = Prob [Ci < Ai < C2, Ci < A2 < C2, . . . , Ci < Aat < C2] = Prob [An,;,, > Ci, A,„ax < C2] . (66) 

In other words, Rn{Cij C2) provides the joint probability distribution of the minimum and the maximum eigenvalue. 
By definition, 

i?Jv(Ci, C2) = f^' ■■■ r -^(^1' ^2, ■ • • , Ajv) dAi dX2... d\N (67) 
"'Ci "'Ci 

where P is the joint pdf in Eq. p2p . As in the previous section, this multiple integral can be written as a ratio of 
two partition functions 

i?iv(Cl,C2) = ^^4^ 

where 

' N 



^n{CiX2) 



K2 K2 




/ . . . / exp 









^A?-^ln(|A.-A, 



dXi d\2 . . . cLXn ■ (69) 



and Zjv(— 00) = r27v(^oo,oo) is the same normalization constant as in Eq. (llSp . Thus, ^n{Cit^2) represents the 
partition function of the Coulomb gas that is sandwiched in the region [Ci,C2] bounded by the two hard walls at its 
boundaries. 

We then evaluate this partition function in the large N limit using the saddle point method. The formalism is 
exactly same as in the previous section. We first define a counting function p{^) that is nonzero only in the region 
< M ^ ^iid is zero outside. The rest of the calculation is similar as in the previous section, except that all 



the integrals run over the region ^ G [21,^2] where zi = Ci/vTV and Z2 = C2/vN. The action T,[p] is exactly as in 
Eq. psp . Thus the partition function, in the large N limit, behaves as 

(Ci , C2 ) = exp S [pj +0{N)] (70) 



where the saddle point density Pc{^^), in terms of scaled variables zi = (i/VN and Z2 — C2/VN, satisfies the integral 
equation 

P^vTdp'^, (71) 

where V indicates the Cauchy principle part. Next we introduce the shift 

p = zi+x (72) 

and define W ^ Z2 — z\. Since z\ < p. < Z2, it follows that < x < W . Thus x — denotes the location of the left 
barrier and x — W denotes the location of the right barrier. In terms of the shifted variable, we rewrite the density 
field as 

Pc{p = zi+x) =^ f{x;zi,W). (73) 
Eq. (|7T|) then reduces to the integral equation 

x + z.^Vrdx' (74) 

JO X — x 

The integral equation ([71)1 can again be solved using Tricomi's theorem. Note that this equation has almost similar 
form as Eq. (|44p except that the integral on the rhs of Eq. (ffH) runs up to W. From the solution of Eq. (H^ 
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presented in Eq. (|47)) we learned that the density / is nonzero only for < a; < L{z) and is zero for x > L(z) where 
L{z) = I [Vz^ + 6 — z] . So, comparing to Eq. ([74|) we see that there are two possibilities: 

(i) li W — Z2 — zi > L{zi), the solution f{x;zi,W) will be exactly the same as f{x;zi) presented in Eq. (|47p . 
In this case, the Coulomb gas does not feel the presence of the right barrier at Z2- In other words, the solution 
f{x;zi,W) = f{x;zi,oo) is completely independent of W and one can effectively put W — > oo, i.e., put the right 
barrier at infinity. Thus in the case, the charge density diverges as x^^^'^ at the left barrier and vanishes at x — L[zi). 

(ii) li W — Z2 ~ zi < L{zi), then the solution will be given by Eq. (|46|) with L replaced by W. Using g{x) — x + zi 
in Eq. (|46p and performing the integral on the rhs we get 



f{x;zi,W)^ 



8TTy^x{W~x) 



[W^ + 4:W{x + zi) - 8x{x + zi) + B'] 



(75) 



where B' is an arbitrary constant. The normalization condition, f{x; zi, W) dx — I, fixes the constant B' — 8. In 
this case, the charge density diverges (with a square root singularity) at the locations of both the left barrier {x = 0) 
and the right barrier {x = W). 

Thus, putting (i) and (ii) together, we find that the solution for the equilibrium charge density /(x; zi, W) for the 
Coulomb gas sandwiched between two barriers is given by 



fix;zuW) 



Stt-^/ x{l — x) 



[l"^ + Alix + zi) - 8x{x + zi) + : 



for < a; < 



(76) 



where 



mm 



W^Z2 



zi, L{zi) = - 



z? + 6 - zi 



(77) 



Clearly, in the limit W oo, we indeed recover the results of the previous section. Summarizing, if one fixes the 
left barrier at zi and varies the position of the right barrier Z2 (equivalently by varying the distance W = Z2 — zi 
between the two walls), one finds that the charge density at the left barrier always diverges. On the other hand, 
the behavior of the density near the right barrier undergoes a sudden change as W increases beyond a critical value 

zi^. The density at the right wall diverges as long as < Wc, i.e., Z2 < zi + L{zi). 

L{zi), the charge density goes to zero at the right edge of the support 



zf + 6- 



But when W > Wc or equivalently Z2 > zi 



at L(zi) < W. 

Having determined the charge density Pc{^j) = f{x = fi+zi] zi, W), the saddle point action S[pc], is then determined 
via the following equation that is analogous to Eq. (|55|) 



npc] 



1 



1 



dx \n{x) f{x; zi,W) + - / dx {x + zi) f{x;zi,W) 



(78) 



where W = Z2 — z\ and /(x;zi,iy) and I are given respectively in Eqs. ((76)) and ((77)) . Denoting the saddle point 
action S(z\, = S[/9c] and evaluating explicitly the integrals in Eq. ((78)) we get for W < L{zi) 



32 ln(2) - 16 In(VK) + 16 zf + GW^ + 16 Wzi 



-2W^zf 



^ 2 W^zi ~ —W^ 
16 



(79) 



For W > i(zi), S{zi,W) becomes independent of W and sticks to its value ^(zi, L(zi)). On the other hand, we 
know that when W ~> oo, 5(zi,oo) must be equal to the action S{zi) for a single wall as given in Eq. ((56p . Indeed, 
one can check explicitly that S{zi, L{zi)) = S'(zi), thus confirming the expectation. 

The partition function, in terms of the scaled variables zi and Z2, then follows from Eq. (|70p 



nN{zi,Z2)^exp [-PN^ S{zi,W) + 0{N)] . 



(80) 



Note that the denominator Zn{—oo) in Eq. ((^7)) is still given by Eq. ([57]) where S{-V2) = (3 + ln(4))/8. Hence 
taking the ratio in Eq. ((HTj) and using Eq. ((501) we get the joint probabihty Rn{Ci,C2) for large N 



Rn{CiX2) = exp 



-PN^-q, f^,^) +OiN) 



(81) 




FIG. 3: The charge density f{x;zi,W) for zi — and W = 1 and W = 2 respectively. When the left wall location zi = 0, 
-L(O) — \/8/3 — 1.63299. Thus the critical value of the distance W between the walls is Wc = 1.63299. On the left panel, 
W = 1 < Wc (subcritical) where the charge density diverges (square root singularity) at the location of second wall x = W. 
On the right panel, W = 2 > Wc (supercritical) where the charge density goes to zero as a; — > Wc = L{0) < W . 



where 

^^,{z,,z,)^Siz,,W = Z2-z^)-^±^ (82) 

with S{zi,W) given by Eq. (j79p . One can check easily that when the second wall moves to infinity, i.e., Z2 oo, 
vE'(zi,oo) = 0{zi) where 9{z) is given in Eq. ([55]) . Thus, Eq. ([5T|l for the joint distribution of Amin and A,nax is a 
generalization of Eq. (|58p that describes only the distribution of Amin ■ 



IV. NUMERICAL RESULTS 



The reader will realize that the numerical confirmation of the analytical results of the previous section is a delicate 
and potentially computation intensive task. For simplicity we will restrict our selves to the case of the GOE but 
the methods used can be extended to the other ensembles. The simplest way to compute the probability that all 
eigenvalues are greater than some value is to numerically generate matrices from the required ensemble, diagonalize 
them and then count the number m+ that satisfy the eigenvalue constraint required. However because of the order 
suppression of this probability found here, for large N the number of matrices m that one would need to generate 
before seeing a single matrix satisfying the constraint is huge. The estimate for the probability that all eigenvalues 
are positive in this method is given by 

Qjv(O) = (83) 
m 

In [3] an approximate argument for the GOE (/? = 1) was made yielding 9—1/4: and a subse quen t numerical study on 
matrices up to 7 x 7 with an A^-dependent fit = aN" yielded a = 2.00387 and a — 0.3291 [2^ was found. However 
given that there are 0{N) corrections and the size of the systems studied are so small this fit cannot be taken too 
seriously. In fact one can use the Coulomb gas representation of the eigenvalues of Gaussian ensembles in order to 
numerically compute 9(0) for much larger values of N. However as a test of this method for smaller values of N we 
may adopt the direct enumeration approach of Q but slightly improve it to gain a few extra values of N . 
If the all the eigenvalues of a matrix are positive then for any vector v we must have that 

(v, Mv) > (84) 

In particular if we choose the vector v to be one of the N basis vectors then Eq. (jSl]) implies that 



(ei,Me,) = Mu > 0, 



(85) 
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and so if M is positive (in the operator sense), all of the diagonal elements must be positive. The estimation of Q7v(0) 
can thus be slightly improved by increasing the chances of seeing a positive matrix by forcing the diagonal elements 
to be positive. With respect the the simplest form of enumeration, the matrices M* generated are the same as those 
for the GOE but the diagonal elements are replaced with their absolute value. The probability that a given matrix 
M has all diagonal elements positive is 1/2^, thus if m*^ denotes the number of these so generated matrices (with 
positive diagonal elements) then the estimation of the probability that a GOE matrix is positive is given by: 

Qn{^)^^ — - (86) 



2^ 



TO 



For small N this method thus appreciably increases the probability of generating positive matrices and thus enhances 
the accuracy of the estimate for (5jv(0). Even so using this method it is virtually impossible to obtain meaningful 
results for N > 8. The results obtained by this modified enumeration method are shown on Fig. ^ (squares). 

For large N it is in fact much better to evaluate Qn{0) directly from Eq. (jlSp via a Monte Carlo method. We note 
that we can write 

^ Z(-cx3) = (G(A)), (87) 



(27r) 

where G(A) = nj:<j 1-^* ~ ^^"^ angled bracket indicates the average is over A.^ taken to be independent and 
Gaussian of zero mean and unit variance. The term Z{0) can also be related to the expectation over A^ which are 
similarly independent and Gaussian of zero mean and unit variance but conditioned to be positive, and we denote the 
average of with respect to these variables by ( • ) + . We find that 

Z{0) = (G(A))_ 



(27r) 

the left-hand side now has has the form of a conditional average, the factor of 2^ giving the correct normalization for 
the conditioned probability distribution. Putting all this together yields 

1 (G(A)) + 

The two expectation values can be computed via Monte Carlo sampling, the unconditioned one by using Gaussian 
random variables and the second trivially by using the absolute value of Gaussian random variables. Shown in Fig. 
(HI) is a plot of ln((5Ar(0)) for the GOE (/3 — 1) ensemble measured as described above (circles), we see that the 
agreement for small N with the results obtained by modified enumeration approach is excellent. For larger values of 
N we have used 5 x 10^ Monte Carlo samplings and there is a significant amount of fluctuation as indicated by the 
error bars. The Monte Carlo results were fitted using the fit a fit ax^ +bx + c, for values of N between 3 and 35, the 
fit yields a = —0.272, b = —0.493 and c = 0.244 which is in good agreement with that predicted here. However given 
the errors for large iV and the fact that there are probably corrections of 0(ln(A^)), the numerical estimate for the 
exponent probably on has about a 10 % accuracy. 

Also, by direct sampling over Gaussian matrices, one can numerically evaluate the the rescaled density for states 
for matrices having only positive eigenvalues. Because we use the direct sampling method we are clearly restricted to 
small values of iV, however in Fig. ^ we show the analytical large N result for / with that computed numerically 
for matrices with N = 6 and for N = 7, we see that despite the small value of N the agreement is already rather 
good, the main deviation being in the tails for large /z. 



V. CONCLUSIONS 



In this paper we have shown how the Coulomb gas formulation of the distribution of eigenvalues of {N x N) Gaussian 
random matrices can be exploited to derive exact asymptotic results concerning the extreme value statistics of their 
eigenvalues. Our main results are summarized as follows. 

(i) We have shown the probability P/v that all eigenvalues are positive (or negative) (or equivalently the probability 
that Amin > or A^ax < 0) decays as Pn ~ exp[-/36'(0) N"^] for large where 6(0) = ln(3)/4 = 0.274653 ... and /3 
is the Dyson index. 

(m) More generally, we have computed the probability Prob[Amax £ t, N] that the maximal eigenvalue is located deep 
within the Wigner sea region, far to the left from its average value V2N, i.e., when t ~ 0{N^^'^) < \/2N. This 
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_400 

10 20 30 40 

N 

FIG. 4: Monte Carlo computation of ln(Qjv(0)) for the GOE (black circles error at large indicated by the size of the circles) 
along with quadratic fit (solid line). Shown as squares are the results obtained by modified enumeration. 



probability has the asymptotic form, ~ exp — /3iV^ <i> (^ ^^^^ * ^ for large N, where the large deviation function <&(?;) 
has been computed exactly. 

(iii) We have also computed the asymptotic joint probability distribution of Xmin and Amax- 

Our result in (ii), valid when ^/2N — t ~ 0{N^^'^) (deep inside the Wigner sea) is complimentary to the Tracy- 
Widom result that concerns the distribution of Amax about its mean value (near the edge) over a small range of 
width ~ N~^/^, i.e., for \/2N — t ~ 0{N~^^^). We have demonstrated explicitly how these two results match up 
smoothly as one approaches from deep inside the Wigner sea to its right edge. 

The key step in our method for computing the distribution of Amin (or Amax) in (ii) consists in using a functional 
integral approach to study the Coulomb gas representation of the problem and imposing a single hard wall constraint 
which enforces the fact that no eigenvalues can be to the left (or right) of a given point. For the computation of the 
joint distribution of Amin and Amax in (iii), we needed to confine the Coulomb gas within two hard walls. In the limit 
of large N , the functional integrals can be evaluated by the saddle point method and the resulting integral equations 
for the saddle point density can be solved explicitly using Tricomi's theorem (23| . 

Our method is actually rather general and has already been adapted to study the critical points of Gaussian 
random fields in large dimensional spaces (lol . [2^ and the extreme value statistics of the maximum eigenvalue of 
Wishart random matrices [2l|. One can possibly find further applications in related statistical problems. For instance 
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FIG. 5: The analytical large A'^ formula for f{x;0) with z — (solid line) along with the numerically generated averaged 
histogram of 6 x 6 (open squares) and 7x7 (solid circles) Gaussian matrices with positive eigenvalues. The agreement is 
already good, the main difference occurring at the large x tail. 



the method is probably adaptable to study the statistics of the index (the number of negative eigenvalues) 6] of 
random matrices and one could also study the extreme value statistics of the minimal value of the modulus of the 
eigenvalues by introducing the the appropriate constraint on the density of eigenvalues in the functional integral 
formulation of the problem. 

Note that in this paper we were able to compute only the leading large N behavior of the distribution of extreme 
eigenvalues. It would be interesting to compute the sub-leading corrections to this leading behavior. Some recent 
attempts have been made in this direction [26|. 
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